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We study the ground state and low energy excitations of vortices pinned to columnar defects in 
superconductors, taking into account the long-range interaction between the fluxons. We consider 
the "underfilled" situation in the Bose glass phase, where each flux line is attached to one of the 
defects, while some pins remain unoccupied. By exploiting an analogy with disordered semiconduc- 
tors, we calculate the spatial configurations in the ground state, as well as the distribution of pinning 
energies, using a zero-temperature Monte Carlo algorithm minimizing the total energy with respect 
to all possible one-vortex transfers. Intervortex repulsion leads to strong correlations whenever the 
London penetration depth exceeds the fluxon spacing. A pronounced peak appears in the static 
structure factor S(q) for low filling fractions / < 0.3. Interactions lead to a broad Coulomb gap in 
the distribution of pinning energies g(e) near the chemical potential (i, separating the occupied and 
empty pins. The vanishing of g(e) at leads to a considerable reduction of variable-range hopping 
vortex transport by correlated flux line pinning. 



PACS numbers: 74.60. Ge, 05.60.+W 



I. INTRODUCTION 

Recently, the static and dynamic properties of flux 
lines in high-temperature superconductors subject to an 
external magnetic field B have been intensively studied 
both experimentally and theoretically Q. From a the- 
oretical point of view, the importance of thermal fluc- 
tuations H and the strong influence of point defects 
pose a variety of interesting problems, leading to 
strikingly rich and complex phase diagrams j^j. For the 
purpose of applying (high-temperature) superconductors 
in external magnetic fields, an effective vortex pinning 
mechanism is essential, in order to minimize dissipative 
losses caused by the Lorentz-force induced movement 
of flux lines across the sample. Therefore, in addition 
to point defects, the influence of extended or correlated 
disorder, promising stronger pinning effects, on vortex 
transport properties has been considered. Experimen- 
tally, it has been found that by high-energy ion irradi- 
ation linear damage tracks are formed in the material. 
These columnar pins bind the flux lines very strongly, 
thus significantly increasing the critical current and shift- 
ing the irreversibility line upwards On the theoret- 
ical side, the pinning of flux lines to linear defects and 
the ensuing transport properties have been studied by 
Lyuksyutov, @] and in detail by Nelson and Vinokur 
In Ref. |^] , correlated pinning by twin and grain bound- 
aries is discussed as well, a topic which has been more 
thoroughly investigated recently by Marchetti and Vi- 
nokur H. 

Similar to the physics of flux lines in pure systems, 
|^| the statistical mechanics of vortices interacting with 
columnar pinning centers, which are aligned parallel to 
the magnetic field (Fig. ??), may be mapped onto the 



quan tum mechanics of bosons in two dimensions Q (see 
Sec. HA). At high temperatures, one finds an entan- 
gled liquid of dclocalized vortices, separated by a sharp 
second-order phase transition (corresponding to a bo- 
son localization transition [fl0||) from a low-temperature 
phase of lines strongly attached to the columnar pins. 
This Bose glass phase is characterized by an infinite tilt 
modulus, and was shown to be stable over a certain range 
of tipping angles of B away from the direction parallel 
to the linear pinning centers (transverse Meissner effect) 
II . In addition, the theory also suggests a Mott insulator 
phase at low temperatures, with both tilt and compres- 
sional moduli acquiring infinite values, occurring when 
the vortex density matches the density of columnar pins 
H . The ensuing phase diagram is sketched schematically 
in Fig. ??. Since the Mott insulator phase is predicted to 
occur at low temperatures deep inside a regime with very 
slow dynamics, it may actually be difficult to access for 
kinetic reasons. It is also possible that the Mott insulator 
vanishes entirely for sufficiently strong disorder. Some of 
the predictions of Ref. Q have now been tested and con- 
firmed by experiments on a number of different samples 
pd]-^6[. E.g., the scaling behavior near the boson local- 
ization transition was studied, [ fti"| and in addition even 
some properties of the Bo se g lass itself have been inves- 
tigated to some detail @Hlf| . 

In the Bose glass phase, the linear resistivity vanishes 
for low external currents / « J Cl and the most impor- 
tant mechanism for vortex transport is "tunneling" be- 
tween different columnar defect sites via the formation of 
a pair of "superkinks" i.e.: a fluxon forms a tongue- 
like double kink extending from one pin to another, usu- 
ally a distant defect with nearly the same energy, such 
that the tunneling probability between the columnar de- 
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fects is optimized (Fig. ??). This is very closely related 
to variable-range hopping transport of charge carriers in 
disordered semiconductors, jl7| and leads to the highly 
nonlinear expression [p.8| 



£ = p Jexp[-{E K /T)(J /J)P] 



(1.1) 



for the current-induced electric field £{J) in the limit of 
very low currents J <§C J c (Fig. ??). In Eq. (1.1), po de- 
notes the normal-state resistivity, Ek is a typ ical kink 
ener gy, an d J sets the current scale, see Eq. (|l.2j ) and 
Sec. HID. By assuming a short-range repulsive inter- 



action between the flux lines, Nelson and Vinokur could 
identify p with the Mott variable-range-hopping expo- 
nent in two dimensions, p = 1/3. 

Similarly, for the case of parallel planar defects con- 
taining the B axis, Marchetti and Vinokur have found 
an analagous formula, with p = 1/2, demonstrating that 
twin and grain boundaries serve as even better pinning 
centers for flux lines . Returning to linear defects, Hwa, 
Le Doussal, Nelson, and Vinokur have suggested that an 
even more drastic reduction of vortex motion may be ob- 
tained by introducing a controlled splay, i.e., dispersion 
in the orientation of the columnar pins |l9| . A new low- 
temperature phase, the "splayed glass" , is found which is 
related to the Bose glass, but has a finite equilibrium tilt 
modulus (however, the system is dynamically frozen due 
to diverging energy barriers preventing the relaxation of 
small externally induced tilts). Flux-line transport in 
this phase should be more strongly inhibited than in the 
Bose glass for two reasons: (i) variable-range hopping 
must now proceed via a much slower process of selecting 
both energy and angle; as a consequence, the ex pon ent 
in the nonlinear current-voltage characteristics jl.lj ) is 
found to be enhanced to p — 3/5 [jl{|. (ii) Flux lines 
in the splayed glass ground state are entangled and can 
only move by cutting through each other, at the expense 
of a considerable amount of (free) energy |2(|j . 

In the above studies, j^j9,[l{| the intervortex repul- 
sion was only treated approximatively, using order-of- 
magnitude estimates. Specifically, it was assumed as be- 
ing essentially short-range. However, in fact the inter- 
action range is set by the London penetration depth A, 
which diverges at the critical temperature, and thus the 
intervortex repulsion may effectively become very long- 
range, namely when A is large compared to the average 
spacing between flux lines, A > a = (4/3) ^(^/-B) 1/2 
(here (f>o = hc/2e denotes the magnetic flux quantum). 
These strong long-range interactions will then determine 
the width and shape of the distribution of pinning ener- 
gies, which, as has been emphasized by Gurevich, pi I in 
turn strongly influences vortex transport properties! Bfl . 

Simultaneous observations of both columnar defects 
and magnetic flux lines, using scanneling tunneling mi- 
croscopy, have recently been reported for NbSe2 sam- 
ples in the "underfilled" regime, where the number of 
columnar pins Nn exceeds the number of vortices N 
|[12| . For the high-temperature superconducting material 



Bi2Sr2CaCu208+5 (BSCCO), a combined chemical etch- 
ing and magnetic decoration technique has been similarly 
successful in determining the positions of both the colum- 
nar defects and vortices; and indeed cases have been 
found, in which the strongly correlated spatial distribu- 
tion of the flux lines suggests that the intervortex repul- 
sion clearly dominates in magnitude over the unavoidable 
fluctuations in pinning energies stemming from the varia- 
tions in the ion track diameters |D| . This would already 
imply a noticeable reduction of vortex transport; for the 
width of the distribution function g(e) of the vortex pin- 
ning energies e, is then no longer determined by the com- 
paratively small variations in pin diameters, but by the 
much larger typical interaction scale. Thus the value of 
g(e) at the chemical potential has to be reduced, leading 
to a greater value for the current scale in Eq. (p~l|), 



J = c/ (j) Q g(p,)d 3 



(1.2) 



where d is the average defect distance. Therefore, by 
a suitable "tailoring" of g(e), the prefactor of the power 
l aw in the exponent of the current -voltage characteristics 
(1.1) may be adjusted in order to achieve more effective 
pinning. 

Moreover, by further exploiting the analogy to two- 
dimensional bosons localized at randomly distributed de- 
fect sites, one has to at least consider the possibility 
that these spatial correlations may produce a "Coulomb" 
gap [fl7f in the distribution of pinning energies g(e) near 
the chemical potential //, which separates the filled and 
empty energy levels. In the limit of infinitely long-range 
interactions, A — > oo, one would expect g(e) to vanish 
near /i according to a power law 
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(1.3) 



see Sec. Ill B| . This of course invalidates Eq. (|l.l[), be- 
cause there a finite value of g{p) ^ has been implicitly 
assumed Instead, a simple re-analysis of the opti- 
mization procedure leading to Eq. (IT) shows that the 
exponent po = 1/3 is to be replaced by the larger value 
p = (s + l)/(s + 3], thus strongly suppressing vortex mo- 
tion (see Sec. [II D| ; also the current scale Jo becomes 
modified). 

In order to determine the gap exponent s, and thus p, 
and also to understand down to which values of X/d and 
ao I d the intervortex repulsion may still be considered as 
sufficiently long-range to produce a correlation-induced 
Coulomb gap, one has to resort to numerical simulations 
in the spirit of previous investigations of the so-called 
Coulomb glass (localized charge carriers interacting with 
a 1/r potential) p2j p4|. Except for a mean-field anal- 
ysis for the Coulomb potential problem, |l~7| the only 
related prior investigation appears to be a scaling anal- 
ysis of logarithmically interacting vortices in two dimen- 
sions by Fisher, Tokuyasu, and Young ]25| ], We have 
thus performed an extensive Monte-Carlo study of a two- 
dimensional system of N localized (i.e., classical) parti- 
cles subject to a Bessel function potential, representing 
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the interacting flux lines, as a function of the two pa- 
rameters X/d and / = N/Nd = d 2 /og. We find that for 
fillings 0.1 < / < 0.4 (d « ao in this range), and large 
values X/d > 5 a very wide and pronounced gap appears 
in the distribution of pinning energies, described by an 
effective gap exponent assuming values up to s c ff ~ 3, 
and as a consequence the transport exponent reaches at 
maximum p e g w 2/3. This means that correlation ef- 
fects due to long-range repulsive forces may drastically 
enhance the correlated pinning of flux lines in the Bose 
glass phase. However, even for the case A « d w ao , one 
still finds a marked minimum in g(e), with an effective 
gap exponent s c s ~ 1, and hence p c g ~ 1/2. Only for 
X/d < 0.5 is the non-interacting result po = 1/3 recov- 
ered. We conclude that the correlation effects due to the 
intervortex repulsion are surprisingly strong, and may 
actually be an interesting means to "tailor" the distribu- 
tion of pinning energies in ion-irradiated samples, such 
that flux creep is suppressed as effectively as possible. 
A preliminary account of some of these results and their 
connection with recent experiments |l3| ] has been given 
in Ref. |(|. 

We remark that we do not expect that additional point 
defects (e.g., oxygen vacancies) could alter our results 
substantially. First, pinning to point defects is subject 
to much stronger thermal renormalization than pinning 
to correlated disorder. Second, for very low currents, 
these point defects might indeed trap the spreading of 
the double-kink configuration considered above; how- 
ever, this becomes effective only on unphysically large 
length scales (in the km range), [|| and thus is irrelevant 
for realistic samples. Finally, no strong pinning effects 
have been reported in superconducting materials above 
the irreversibility line prior to heavy-ion irradiation. 

This paper is organized as follows. In the subsequent 
section we describe how our effective Hamiltonian may 
be derived from the free energy of interacting flux lines 
subject to columnar disorder, introduce the relevant pa- 
rameters and energy scales, and discuss the region in the 
phase diagram where we expect our model to be valid. 
Furthermore, we describe our specific Monte Carlo sim- 
ulation algorithm, also discussing the different possible 
choices of initial configurations for the energy minimiza- 
tion procedure. In Sec. [I], we describe our numerical 



results and give a number of examples for the spatial 
configurations we find, as well as for the shape of the 
ensuing distribution of pinning energies, as we vary the 
filling / and the interaction range X/d. We also discuss 
vortex transport mechanisms in the Bose glass phase, and 
deduce the form of the current -voltage characteristics in 
the variable-range hopping regime from the previously 
determined distribution of pinning energies, identifying 
the effective Mott exponent p c ff- A brief summary and 
discussion of our results concludes this work. 



II. DESCRIPTION OF THE MODEL AND 
MONTE CARLO SIMULATION PROCEDURE 

A. Model equations 

We start by considering the following model free en- 
ergy for N flux lines, defined by their trajectories Yi{z) 
as they traverse the superconducting sample of thickness 
L, interacting with each other and with Ns columnar 
defects, parallel to the magnetic field B which is aligned 
along the z axis (perpendicular to the copper oxide planes 
in the case of the anisotropic high-T c materials) , |y (see 
Fig. ??) 

F „ l{ri } ] = ^ d2 g{|(^)) 2 + Ig Fhi(z)1 

No ~| 

+ ^F D [r 4 (z)-R fe (z)] I . (2.1) 

k=l ) 

The first term describes the elastic line tension, and stems 
from an expansion with respect to small tipping angles 
of the line energy of nearly straight vortices ||. For a 
dense liquid of flux lines, A > ao, the tilt modulus is 
£i [M±/M z )eo ln(ao/£), where the material anisotropy 
is embodied in the effective-mass ratio M±/M z (<C 1 for 
high-temperature superconductors), £ is the (in-plane) 
coherent length (k = A/£ « 100), and e — (</> /47rA) 2 
sets the energy scale. On the other hand, for low fields 
B < 0o/A 2 , i.e.: A < ao, ?i « eo because of the mag- 
netic coupling between the copper oxide planes 0. In 
this situation the vortex lines are therefore considerably 
less flexible than at high fields. 

Furthermore, (z) = \vi(z) — Tj(z)\, and 



V(r) = 2e Q K Q (r/A) 



(2.2) 



represents the repulsive interaction potential between the 
lines; the modified Bessel function Ko(r/X) describes a 
screened logarithmic interaction, for Kq(x) oc — lnx for 
x — ► 0, and Kq{x) oc x" x l 2 e~ x for x — > oo. Thus the 
(in-plane) London penetration depth A defines the in- 
teraction range. Note that according to the two-fluid 
model, the penetration depth diverges at the zero-field 
transition temperature T c , 



X(T) = Ao [1 - (T/T c f 



1-1/2 



(2.3) 



and therefore the interaction range becomes considerably 
longer upon approaching T c ; e.g., for T/T c « 0.93 one 
has A(T)/Ao = 2. Because eo oc A~ 2 , the interactions are 
both weak and long-range near T c . 

Finally, the columnar pins are described by a sum of 
N-q z-independent potential wells Vd , with average spac- 
ing d, centered on the randomly distributed positions 
{R/c}. The damage track diameters, produced by heavy- 
ion irradiation, are typically 2co ~ 100A, with a varia- 
tion induced by the root mean-square dispersion of the 
ion beam of dck/co ~ 15%. The pinning energies Uk are 
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related to the column diameters via the (interpolation) 
formula Hi 



In 



1 



(2.4) 



the variation of which thus induces a certain distribution 
of the pinning energi es P. with a width w = U (SU 2 ) 
determined from Eq. (2.4), 



w 



1 + (V2£/c ) 2 co 



(2.5) 



(Notice that here the repulsive vortex interactions have 
not been accounted for.) Above a certain temperature To, 
however, the effective radius of the normal-conducting 
region at the defects will be given by the coherence 
length, which grows with temperature according to 



£(T)=6,(1-T/T C 



-1/2 



(2.6) 



hence \/2£(T) should be inserted into Eq. (2.4) instead 
of c k , whenever \/2£(T) > c , or T > T , with 



T Q = T c (1 - 2^/c 2 ) 



(2.7) 



Consequently, above To the effective fluctuations of 
the pinning energies P (without interactions) will be 
smoothened out. For Co ~ 50A and £o ~ lOA, e.g., this 
will happen at To R> 0.92T C , and the characteristic pin- 
ning energies will be Uq — {Uk) ~ 0.35e . On the other 
hand, at T = one finds instead U$ ps 1.3eo, and, as- 
suming Sck/ck ~ 0.15 one gets w ps 0.13eo, while for 
T = 0.8T C one has £ ~ 2.24£ , leading to U w 0.63e 
and w ps O.lleQ. 



The thermodynamic properties of the model (2.1) may 
be derived from the grand-canonical partition function 



of the vortex line problem to the quantum mechanics of 
charged bosons is summarized in Table | (from Ref. ||). 
Notice that in this analogy, the real temperature T as- 
sumes the role of Planck's constant Ti in the quantum 
representation, and the boson electric field and current 
density map on the superconducting current J and the 
true electric field £, respectively; hence the roles of con- 
ductivity and resistivity become interchanged. 

In the following, we shall be interested in the low- 
temperature properties of flux lines in the Bose glass 
phase, pinned to columnar defects, with filling fraction 



/ = N/N B = {d/a ) 2 = B/Bj, < 1 



(2.10) 



More specifically, we need B < B*(T), where vortex in- 
teractions become important in determining the localiza- 
tion length, and thus affect the pinning of the flux lines 
to the columnar defects considerably |8j. In the temper- 
ature range of interest here, 



B*(T) 



B = (4/3) 1 /20 o / d 2 
B^c /2^) 2 {1-T/T c ) 



for T < T 

for T < T < Ti 



(2.11) 



Ti denotes the temperature above which the entropy as- 
sociated with thermal fluctuations becomes relevant for 
the vortex pinning, by modifying the localization length 
and thus renormalizing the binding free energies (for 
the temperature dependence of B*(T) for T > T\, see 
Ref. §, App. D). The estimate of Ref. § is 



Ti «T C 



( Co /4g )(ln K /Gi) 1 / 2 
l + (c /4eo)(ln/t/Gi) 1 /2 



(2.12) 



oc » N 

N=0 ' J t=l 



-F N [{r,( Z )}]/T 



(2- 



where /x denotes the chemical potential 



fi = eo In k 



= ^ {H C1 - H) 
4tt 4tt v 1 ' 



(2.9) 



which changes sign at the lower critical field H, 
4>q In k/AttX 2 . As is explain ed i n deta il in Refs 



the statistical mechanics of (2A) and (2.8) can be for- 
mally mapped onto a two-dimensional quantum mechan- 
ical problem by using a transfer matrix approach, the 
physical z direction being interpreted as an imaginary 
time axis. In the ther mod ynamic limit L — > oo, the prop- 
erties of our model (2.1) are determined by the lowest 
energy eigenvalue, i.e., the corresponding ground-state 
wave function, which is symmetric with respect to ex- 
change of flux lines, and thus of bosonic character; fur- 
thermore the effective temperature of this system of inter- 
acting bosons is then to be taken as zero. The mapping 



where Gi = (M Z /2M ± )(T C / H 2 ^) 2 is the so-called 
Ginzburg number, and the thermodynamic critical field 
at zero temperature is H c = v / 2</>o/47rAo^o- For YBCO, 
typically Gi » 0.01, and with k ps 100 we get T x /T c m 
0.96 and B*{T X ) w O.25S ; for the highly anisotropic 
material BSCCO, upon using Gi ~ 0.1 instead, one finds 
Tx/Tc « 0.89 and B*{T X ) w O.69Ti . For T < T u each 
flux line is localized on essentially one columnar defect, 
and thermal fluctuations play a very minor role; thus all 
temperatures in this range may be considered as "low" . 
At Ti, entropic effects become important, and at the even 
higher so-called depinning temperature Tj p the localiza- 
tion length of a single vortex is determined by the inter- 
play of a large number of pins (see Fig. ??). Similarly to 
Ti, Tip may be estimated as || 



T 



dp 



T r 



(coa/4$ )(ln K /Gi) 1 / 2 
l + (c a/4e )(ln K /Gi) 1 /2 



(2.13) 



where a = ln 1/2 [(d/V2£o)(l - T 1 /T c ) 1 /2]. With d 
lOOOA, e.g., using the above numbers one finds Td p /T c 
0.98 for YBCO, and T dp /T c » 0.94 for BSCCO. 
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For T less than the characteristic fluctuation temper- 
ature T\ , we may thus simply consider the classical limit 
of the corresponding boson problem (Ti — > 0); further- 
more, as the vortices are well separated in this regime 
[with B < B*(T)], the Bose statistics become irrelevant. 
Thermal wandering is effectively suppressed, and the flux 
lines wi ll b e essentially straight; hence the tilt energy 
in Eq. (2.1) can be neglected. Therefore, in the boson 
representation we eventually have to consider the time- 
independent problem of N interacting particles located 
at > N available defect positions. Thus we define 
our model by the two-dimensional effective Hamiltonian 



Nu 



N D 



H=-Y,n i n j V{r ij )+Y j n i t i , 

i=£j i=l 

and its grand-canonical counterpart 



No 



(2.14) 



i=l 



(2.15) 



Here i, j = 1, . . . , Nu denote the defect sites, randomly 
distributed on the xy plane; ni = 0, 1 represents the cor- 
responding site occupation number (n^ = 1 if a flux line is 

bound to columnar defect i), Y^i=\ n i — N. We have also 
included random site energies tj , or iginating in the vari- 
ation of pin diameters [see Eq. fl2.5|) 1. Their distribution 
P can be chosen to be centered at (t) — 0, with width 
w; this amounts to absorbing the average pinning energy 
Uq = (Uk) (which may include small thermal renormal- 
izations) into the chemical potential /i. Realistically, the 
probability distribution of the site energies would likely 
be Gaussian; Jl3|] however, for simplicity we shall assume 
a flat distribution 



P(ti) = G(w-\ti\)/2i 



(2.16) 



[O(x) denotes the Heaviside step function]. As the total 
shape and width of the distribution of pinning energies 
will turn out to be determined by the interactions, the 
precise form of P(U) is of minor importance. Note that 
P corresponds to a bare densit y of s tates. 

For the interacting system (2.14), we define single- 
particle site energies ti as follows: For filled sites (n, = 1), 
ti is the energy required to remove the particle at site i 
to infinity; for empty sites (n, = 0), correspondingly ti 
is the energy needed to introduce an additional particle 
from infinity to site i. Equivalently, we may state that ti 
constitutes the potential energy on site i, resulting from 
the disorder and the interaction with all other occupied 
sites j, or 



No 



(2.17) 



In thermal equilibrium, the chemical potential /j, sepa- 
rates the occupied and empty states, ti < fi for all i 



with rii = 1, and tj > /i for all j with rij = 0. The 
distribution of pinning energies g(e), now with the inter- 
vortex repulsions taken into account (as opposed to P), 
can be viewed as an interacting single-particle density 
of states and may be obtained from the statistics of the 
energy levels tf, i.e., g{t)dt is defined as the number of 
states per unit area with energies in the interval [e, t+dt\. 
Note that in a system of area L\ this corresponds to a 
normalization of the density of states according to 



g(t)dt = N B /Ll = I /d 2 . (2.18) 



In the following, the terms "distribution of pinning ener- 
gies" and "(single-particle) density of states" will be used 
synonymously. For finite A > ao, the chemical potential 
may readily be estimated as 

^^(X/a a fV(a a )^t f(X/d) 2 \n(f\ 2 /d 2 ) , (2.19) 

as each line interacts with approximately (A/ao) 2 other 
vortic es. Actually, the prefactor of the logarithm in 
( |2.19 ) is independent of A, because to oc 1/A 2 (27). Simi- 



larly, the typical overall width 7 of the energy level dis- 
tribution is 



1 ^V{a Q )^t \n(f\ 2 /d 2 ) 



(2.20) 



and roughly the "bandwidth" for the occupied levels will 
be f=a /j. The relevant energy scales in our problem are 
thus ordered according to /i > fj 3> w. 

Another important quantity is the hopping energy 
Ai^j associated with the transfer of a particle from an 
occupied site i to an empty site j (conserving the to- 
tal particle number N). The simplest way to determine 
this energy cost or gain |17| is to decompose the process 
into two steps, namely removing a particle from site i to 
infinity, thus gaining the energy e,-, and then taking it 
from there back into the system at site j, which costs the 
amount tj — V(ry); the additional contribution V(rij) 
here stems from the fact that after removing the particle 
at site i there were only N — 1 particles left, while tj was 
defined for a TV-particle system, and thus the (fictitious) 
interaction with site i has to be accounted for explicitly. 
Hence one finds 



A; - = 



Vinj) 



(2.21) 



and certainly A,-_>j > for all pairs of sites with m — 1 
and rij = is a necessary condition for the ground state 
configuration; in fact, a sufficient condition for a stable 
ground state may be stated by demanding that the ener- 
gies associated with any m-particle hops all be positive. 
Thus a hierarchical scheme for testing specific configu- 
rations can be constructed, [Q in which one first tests 
a putative ground state for stability against one-particle 
hops, two-particle hops, etc. 

With regard to th e reg ime of ap plicability of the ef- 
fective Hamiltonian (2.14) or (2.15), we note once more 
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that for T > Ti entropic corrections become important, 
and for T > Td p a flux line can certainly be no more 
associated with a single columnar defect, and our sim- 
plified model breaks down. However, the above esti- 
mates demonstrate that thermal fluctuations are in fact 
essentially negligible up to temperatures very near the 
irreversibility line. Similarly, in order that all the flux 
lines stay pinned to the columnar defects, the condition 
B < B*{T) must be fulfilled, which provides us with up- 
per limits for the filling /, above which the intervortex 
repulsions will depin some of the flux lines, and move 
them to interstitial sites between the columnar defects 
(compare Fig. ??). 

Finally, we remark that upon introducing the spin vari- 
ables 



ai = 2rij - 1 



(2.22) 



hence <7i = +1,-1 for = 1,0, respectively, the 
effective Hamiltonian ( [2.15 ) is mapped onto a two- 
dimensional random-site, random-field antiferromag- 
nctic Ising model with long-range exchange interactions, 



E 



with 



V{rij) 
-U 



> 



4^ 



vinj) 



and 



b = 55>-a0-s£v(t 



(2.23) 

(2.24) 
(2.25) 

(2.26) 



Finding the density of states for this system by any ana- 
lytical means beyond mean-field type considerations, fl7} | 
or phenomenological scaling arguments, pj| constitutes 
a very difficult task, and has eluded successful treatment 
to date. Therefore we have to resort to numerical stud- 
ies using a suitab le M onte Carlo algorithm; fortunately, 
the Hamiltonian (2.14) is precisely of the form studied in 
the context of charge carriers localized at random impu- 
rities in doped semiconductors (Coulomb glass problem) , 
[ p7| , p2|j24| and we may utilize the considerable experience 
already gathered in these previous studies. 



B. Simulation algorithm 

Our simulation algorithm closely follows the energy- 
minimization procedure as described in detail by Efros 
and Shklovskii in Ref. p7| , Chap. 14. In a square of linear 
size L±, a random number generator chooses Ad defect 



coordinates (by taking L± = \/Ad, the average defect 
distance d is held fixed); we have performed simulations 
for N D = 100, 200, 400, and 800. Then each site is as- 
signed a random number drawn from the distribution 



( |2.16 ), the width of which we have chosen as w/2eo = 0.1, 
and the interaction potentials V(rij) are calculated ac- 
cording to ( |2.2D for all pairs (ij). Following Ref. p2| , we 
have used periodic boundary conditions, thus closing the 
initial square to a torus and thereby reducing boundary 



effects. Here, 



is defined as the minimum value of the 



distances between site i and the nine equivalent sites j 
in the original and repeated "cells" . For the interaction 
range A, we have studied the values A/rf = 0.5,l,2,5, and 
oo. In order to handle the latter case of infinite-range, 
purely logarithmic interaction carefully, we have applied 
the Ewald sum technique, which amounts to including 
the interaction of each particle with all its periodic im- 
ages; this is achieved by performing the corresponding 
series in part in direct, and in part in Fourier space. For 
details, we refer to Ref. ^8|, App. A. The data ti and 
V{rij) can now be stored for later use. 

In the next step, N = /Ad (/ < 1) of these sites i 
are assigned the occupation number = 1. We have 
employed two very different schemes to generate this ini- 
tial configuration: (i) the occupied sites were chosen ran- 
domly, and (ii) a suitable triangular lattice, with a num- 
ber Ai at of sites as close as possible to the correct A, was 
superimposed on the two-dimensional array of defects, 
and then deformed until the former lattice sites matched 
A^t of the defect positions; the remaining N — N\ at sites 
were then filled at random. The choices (i) and (ii) nat- 
urally correspond to minimal and maximal initial spatial 
correlations, respectively. Using Eq. (2.17), now the site 
energies ej are calculated. Almost certainly, the initial 
configuration will correspond to a highly non-equilibrium 
situation, in the sense that many of the occupied sites i 
will have higher energies than a large number of empty 
places j. In order to relax this configuration by successive 
particle-hole transitions, the occupied site p of highest 
energy, e p = ruax/ ni= i}ej, and the empty site q of lowest 

< e p , the 



energy, e q 



l{n i= 0} e J 



are determined. If e„ 



particle at site p is moved to site q (the corresponding oc- 
cupation numbers are interchanged) and the site energies 
are recalculated. This procedure is repeated until even- 
tually e p < e q ; now the occupied and empty states are 
separated, and the chemical potential can be (approxi- 
mately) obtained from fi = (e p + Cq)/2. 

However, this intermediate state is still a very poor 
approximation to the real ground state, as in general 
it will be unstable towards single-particle hops. In or- 
der to test this, the hopping energies Aj_>j are cal- 
culated for all possible pairs of occupied sites i and 
empty sites j, using Eq. ( [2.21 ). Among those values of 
Aj_>j, which are negative, the minimum is searched for: 
Ap_» 9 = min{Aj_j<o} A,_>j ■. Next the particle transfer 
from site p to site q is performed, thus reducing the to- 
tal energy of the system. Then the new site energies are 
calculated, and in general the "equilibration step" of the 
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previous paragraph will have to be repeated. Once more 
all the A,-_>j are determined, and the whole procedure is 
done again, until finally all the hopping energies acquire 
positive values. The ensuing configuration is accepted as 
a fair approximation to the real ground state for the par- 
ticular distribution of defects, filling fraction /, and inter- 
action range X/d. The algorithm thus yields the positions 
of the N flux lines, the corresponding site energies, and 
the approximative chemical potential fi. By repeating 
this procedure for a number k of different initial config- 
urations (we took k = 16000/iVrj), one may then obtain 
an ensemble-average of the chemical potential (fi) , of the 
density of states (g(e)} (from the site energy statistics), 
and similarly averaged correlation functions, etc. 

Of course, to find the correct ground state, one would 
in principle have to test each configuration against any 
simultaneous m-particle hops, m = 2, ... ,00. However, 
previous investigations have shown that terminating at 
m = 1 already yields at least qualitatively reliable esti- 
mates for the statistics of the energy level distributions 
jl7],^2[ 24| . In fact, the "true" ground state may be very 
difficu t to reach for the corresponding real system as 
well, and the static and dynamic properties may at least 
at low temperatures be satisfactorily described by these 
"pseudo" ground states. The quantitative reliability of 
the results can furthermore be tested by studying vari- 
ous system sizes, and by different sampling procedures 
for the site energies. E.g., one can determine (<?(e)) and 
(y) by directly averaging the energy level statistics of the 
k different realizations, or by performing the average on 
the relative energies e = e — \i calculated in each of the 
realizations separately. It turns out that these two ver- 
sions yield practically identical distributions of pinning 
energies, the tiny differences being confined to the imme- 
diate vicinity of (/x) , and readily interpreted as finite-size 
effects (the second procedure always leads to g(e) = 
for |e| = |e — /i| < V(L±), because the lowest possible 
value of e is of the order V(L±), see Ref. fn]].) One 
can also check that the two very different initial condi- 
tions described above lead to similar final configurations. 
In Sec. HID, wc shall describe how the essentials of the 



current-voltage characteristics in the variable-range hop- 
ping regime may be obtained from the thus determined 
single-particle density of states. (We shall drop the ex- 
plicit notation (. . .) for ensemble averages from now on.) 



III. STATIC CORRELATIONS AND TRANSPORT 
IN THE BOSE GLASS PHASE 



A. Spatial correlations 

We begin the discussion of our results with a descrip- 
tion of the spatial distribution of flux lines among the 
underlying randomly distributed columnar pins. Fig. ?? 
presents the final configuration in the xy plane (perpen- 
dicular to both the magnetic field and the defect ori- 



entation) for three typical simulation runs, performed 
on the identical defect configuration (N^ = 400), with 
infinite-range logarithmic interaction (A — > 00), a nd us - 
ing site energies drawn from the flat distribution (|2 . 1€ ) 
with w/2eo = 0.1 (compare the estimates in Sec. HA), 
but with different filling fractions / = 0.1, / = 0.2, and 
/ = 0.4, respectively, in each case starting from a ran- 
dom initial distribution of vortices (the final states stem- 
ming from the other extreme of possible initial condi- 
tions, namely a distorted triangular lattice, are not iden- 
tical but look qualitatively very similar, see below). In 
contrast to the Poisson distribution of defects (open cir- 
cles), showing the characteristic voids as well as cluster- 
ings, the flux lines (full circles) form a highly correlated 
structure, trying to keep as distant as possible from each 
other to minimize the interaction energy. For low fillings 
[Figs. ??(a) and (b)] the ensuing spatial distribution of 
flux lines resembles a distorted triangular lattice. How- 
ever, for higher filling fractions, / > 0.4 say, the vortices 
are increasingly forced to accomodate the underlying ran- 
dom pin distribution. As can be seen in Fig. ??(c), they 
then tend to occupy the energetically favorable sites en- 
circling spatial voids in the defect distribution. Upon in- 
creasing / even more, beyond / > 0.7 say, (which may be 
achieved by applying higher magnetic fields B) the spa- 
tial correlations disappear on longer length scales, and 
the only remaining effect of the intervortex repulsion is 
to prevent simultaneous occupation of nearest-neighbor 
sites. One has to keep in mind, however, that in this 
regime our basic assumption that each flux line be at- 
t ache d to a single defect, which requires B < B*(T) 
[Eq. ( 2.11 )1, is be ginn ing to break down. As has been 
explained in Sec. II A, for B > B*(T) interactions will 
become important in determining the vortex localization 
length, and eventually force some of the flux lines to leave 
the columnar defects and move to interstitial sites. Most 
of our numerical work is restricted to / < 0.4 for this 
reason. We note that all of these features do not very 
strongly depend on the actual value of the parameter 
X/d, as long as X/d > 1, which ensures that the typical 
interaction energies exceed the fluctuations in pinning 
potential depths w. 

More quantitative conclusions can be drawn from the 
calculation of the static structure factor 5 f (q), 



1 N 



= iv> (q) 



(3.1) 



which is readily obtained from the two-dimensional 
Fourier transform n(q) = e iqFi f the vortex density 
n(r) = J2i S(r — ri); more precisely, we shall consider the 
average of (3.1) over the directions in momentum space, 



S(q) 



dn^ S(q) 



271 



d<t>S{q,4>) 



(3.2) 



In Fig. ??, the calculated S(q) is depicted for both the 
initial and final vortex configurations, for filling fractions 



7 



/ = 0.1 (a), and / = 0.2 (b). These curves were ob- 
tained from averaging over k — 40 different runs with 
Ns = 400, w/2eo — 0.1, and X/d = 5. Again, the corre- 
sponding pictures for any X/d > 1 (keeping w fixed) look 
very much alike. Note that the peak in the dashed curve 
in Fig. ??(a) at qao « 47r/\/3) i.e., at the Bragg peak 
corresponding to the triangular lattice, is shifted to lower 
values and diminished in height in Fig. ??(b); for higher 
fillings the underlying random pin distribution enforces 
an average separation of vortices larger than the ideal 
hexagonal close-packed geometry would permit. For sim- 
ilar reasons, the peak of S(q) for the equilibrated config- 
urations (corresponding to occupied sites in Fig. ??) is 
displaced to lower values with respect to the ideal tri- 
angular lattice, namely to qao ~ 0.86(47r/v3), and also 
broadened considerably due to the appearance of a large 
number of topological defects (see Fig. ?? below). In- 
terestingly, for / = 0.1 there is a double-peak structure. 
The overall appearance, width and peak height of the 
dotted and solid curves in Fig. ?? are very similar, sug- 
gesting the relative unimportance of the specific choice 
of initial conditions. In addition to the suppression of 
fluctuations at low qao < 4 and the marked peak in S (q) 
at (roughly) qao sa 2ir, the shallow dip in the interval 
9 < qao < 12 should be noted. For these spatial cor- 
relation effects, the most sensitive parameter turns out 
to be the filling fraction /, as illustrated in Fig. ?? for 
A/d = 2. While for / = 0.1 and / = 0.2 the structure 
factor is practically indistinguishable from the previous 
pictures for A/d = 5, already for / = 0.4 the peak height 
becomes considerably reduced, and the shallow dip has 
disappeared. For / = 0.8 the peak has vanished, and 
even the correlation gap at low q is halfway closed, thus 
rendering S(q) more and more uniform as for a Poisson 
distribution (cf. the long-dashed curves in Fig. ??). 

Returning to the configuration plots in direct space 
(Fig. ??), we can analyze the final structures more thor- 
oughly by performing a Delaunay triangulation, thereby 
exposing the topological defects (with respect to an ideal 
six-fold coordinated triangular lattice) , simply by deter- 
mining the ensuing coordination number iV c at each vor- 
tex site. The results of this procedure for the configu- 
rations of Figs. ??(a)-(c) are shown in Fig. ??. Clearly 
the resulting vortex configuration can only be viewed as 
a highly distorted "triangular lattice", with its orienta- 
tional correlations being strongly diminished. 

In addition to these gross interaction effects, there 
emerge subtle correlations in the single-particle energies, 
namely the spatial clustering of those sites whose ener- 
gies ti differ very little from the chemical potential /i. 
This remarkable effect has previously been observed by 
Davies, Lee, and Rice for the Coulomb glass case j22| . Wc 
illustrate this property for the distribution in Fig. ??(c), 
i.e.: N D = 400, w/2e = 0.01, A -> oo, and / = 0.4. In 
Fig. ??, along with the entire configuration, the spatial 
distribution of those filled and empty sites with energies 
in the range |e< — fj,\ < 5 are shown, with S/2eo — 0.6 (a), 
and S/2eo = 1.2 (b), respectively [compare Fig. ??(b) 



below]. The occupied sites with energies close to /i are 
strikingly clustered in space, and the corresponding low- 
energy empty sites seem to occupy the complementary 
region. The interactions cause long-range spatial fluc- 
tuations in the pinning energies, which are apparently 
yet not understood in detail. Of course, the very fact 
that a certain occupied site attains a high site energy 
must mean that there are other occupied sites nearby. 
It is less clear, however, why the energies of nearby oc- 
cupied sites for S -C /i should be so similar, or why the 
states in the low-energy part of the empty pseudo-band 
should be similarly correlated in space. Note that the 
Bose glass is not symmetric with respect to exchange of 
particles and holes, since the empty sites do not interact. 
Thus, defining = ej — fj,, those states i,j with both 
and \ej\ < <5, are situated nearby when e^e,- > 0, and are 
widely separated in space when < 0. Upon increas- 
ing the threshold S, both the occupied and the empty 
low— |ej| clusters appear to percolate through the system, 
finally becoming increasingly mixed for 5 > n, and the 
complete configuration is reached. 

Recently, the positions of both columnar defects and 
flux lines have been determined simultaneously in ex- 
periments on NbSe2 |Q and BSCCO samples 13 1, us- 
ing scanning tunneling microscopy and a combination 
of chemical etching and magnetic decoration techniques. 
Indeed, the manifest spatial correlations in Fig. 4 of 
Ref. |b|] strongly suggest the relevance of vortex interac- 
tions. In a recent note, p6[ we used such experimental 
data on the distribution of pins and vortices to infer the 
density of states and the ensuing transport characteris- 
tics, which are thus subject to direct verification by mea- 
surements (the relevant parameters in that case turned 
out to be / w 0.24, X/d « 0.96, and w « 0.1e ). This 
possibility of quantitative comparison of the spatial con- 
figurations and correlations as obtained from simulations 
with those measured in actual experiments opens a fasci- 
nating new field for detailed studies, which have been far 
from feasible in the otherwise closely related disordered 
semiconductor systems [DJ. Even the striking cluster- 
ing in Fig. ?? may possibly be observed experimentally, 
if the current techniques of manipulating individual vor- 
tices with the tip of a magnetic force microscope pjl are 
further refined, which might eventually permit the mea- 
surements of individual pinning energies ej. 



B. Density of states and Coulomb gap 

Before we proceed to present our simulation results for 
the single-particle density of states (i.e.: the distribution 
of vortex pinning energies, taking their interactions into 
account), we briefly summarize Efros and Shklovskii's 
qualitative arguments that long-range repulsive forces 
lead to a soft gap in g(e) [17]]. Let us consider parti- 
cles in D dimensions, which are localized on randomly 
distributed sites r^, and interact via an arbitrary poten- 
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tial 



V{rij) = 2e K( rij /X) , 



(3.3) 



where, quite generally, eo sets the energy scale, A de- 
fines a characteristic length, and K(x) is a positive, 
monotonically decreasing function approaching zero for 
x — > oo. Its inverse therefore exists, and shall be denoted 
by K(x) = K-^x). 

The simplest argument to determine the interacting 
density of states now is as follows (see Chap. 10 of 
Ref. Ml). Consider site energies ej located in an 
interval of width e around the chemical potential fi, but 
situated on opposite sides of [i, e.g. with m = 1 and 
rij = 0. Then in the ground state a particle trans- 
fer from site i to j necessarily cost s the positive energy 

— >j Cj 

V(r l3 ) < Ej 

occur over a distance 



Q - V(rij) > [Eq. ( |2.21fl ], which implies 
-£,*<?. Therefore, hops typically have to 



r y (e)> fl(e)wXA-(e/2eo) 



(3.4) 



The concentration of relevant, i.e.: energetically avail- 
able, impurities hence is n(e) oc R(e)~ D . Thus for small 
e the density of states g(e) oc dn(e)/de becomes 



CD 



\K\x)\ 



2e X D K(x) 



D+l 



(3.5) 



e/2e 



where cd is a proportionality factor of order 0(1). 

We can now specialize to the true long-range poten- 
tials K (x) = x~ a with < a < D and K(x) = -Inx 
(which formally results from the previous case in the 
limit a — > 0). In these situations, the prefactors cd 
can actually be determined using the somewhat refined 
self-consistency argument by Efros pp| , wh ich yields 
c D {a) = 2D(1 - a/D)T{D/2 + l)it- D /\ Eq. @ then 
gives for the inverse-power interactions 



<?(?) 



pr(D/2 + i) fj_ 

air D / 2 X D e \2e 



D/a-l 



(3.6) 



i.e.: the density of states near th e ch emical potential van- 
ishes as a power law [cf. Eq. (1.3)1 with gap exponent 
sq(D) = D/a—1; e.g.: for the Coulomb potential (a = 1) 
in D = 2 and D = 3 dimensions one has sq(2) — 1 and 
so (3) = 2, respectively. For the purely logarithmic inter- 
action, one similarly finds 



9(e) 



DT(D/2 + l) 
w D / 2 X D e 



■ exp 



Di 
2T 



(3.7) 



i.e.: although g(e) does not vanish for e = [/,, there is an 
exponential dip near the chemical potential. Naturally, 
there exists a certain upper cut-off for large e, beyond 
which the density of states assumes its usual shape. 

It is important to note that Eqs. fl3.5|)-(3.7), need not 
be quantitatively correct. Correlations in the distribution 



of site energies have been neglected, and the long-range 
fluctuations embodied in the remarkable clustering de- 
scribed at the end of the previous subsection (Fig. ??) 
are not accounted for. Also, only single-particle trans- 
fers were considered in the above discussion, and multi- 
particle hops could possibly modify the results, especially 
in the case of purely logarithmic interactions pq| . In fact, 
Monte Carlo simulations using variants of the above- 
described algorithm 22 24J yield gap exponents some- 
what larger than the mean-field estimates sq(D). E.g., 
the most detailed study of the Coulomb glass performed 
by Mobius, Richter, and Drittler found s(2) = 1.2 ± 0.1 
and s(3) = 2.6 ± 0.2 in D — 2 and D = 3 dimensions, 
respectively, for the best fits to their results on very large 
systems (up to 40 000 and 125 000 sites, respectively). 
These authors, moreover, were ultimately unable to de- 
cide whether the density of states near the chemical po- 
tential is cor rect ly described by a power law at all. Al- 
though Eq. (3.7) for logarithmic repulsions may have to 
be altered as well, a power-law fit in the energy-range 
accessible to our simulations appears not unreasonable. 

In Fig. ??, we present some of the simulation results for 
the density of states of our effectively two-dimensional 
system with the interaction ( p.2[ ), in the limit A — > oo, 
i.e.: for a purely logarithmic potential. They were ob- 
tained by sampling the site energies of k = 40 runs each 
with either random initial conditions (open circles) or 
starting out with a distorted triangular lattice (filled tri- 
angles) , using a system size of TVd = 400 defect sites and 
a random site energy distribution of width w/2eo = 0.1. 
In this case of infinite-range interactions, we have used 
the Ewald summation procedure in order to obtain reli- 
able results, and before taking the limit A — > oo a term 
oc A 2 has been subtracted from the site energies, ]3l| 

e' l = e l -2^f{X/d) 2 . (3.8) 

The specific choice of the initial conditions for the Monte 
Carlo algorithm has hardly any influence. Although there 
is a general tendency for the random initial configura- 
tions to yield lower total energies on average, this slight 
difference lies well in the range of the actual fluctuations 
in /i over the k = 40 runs with different defect distribu- 
tions. At any rate, the emerging soft "Coulomb" gap near 
the chemical potential is very distinct and surprisingly 
broad for both cases / = 0.2 [Fig. ??(a)] and / = 0.4 
[Fig. ??(b)]. The overall width of the site-energy dis- 
tribution is determined by the interactio ns, an d slightly 
increases with / = N/N D = B/B,p [Eq. (|2.20|) 1. As the 
filling fraction grows, however, the width of the Coulomb 
gap itself (at half maximum, say) does not change. Its 
relative width thus decreases sowewhat with larger fill- 
ings, as it "sweeps through" the density of states along 
with the chemical potential. More important, the effec- 
tive gap exponent apparently depends monotonically on 
/, attaining its largest values for small fillings. 

This becomes more clearly visible in the double- 
logarithmic plots of Fig. ??, where both the filled (full 
symbols) and empty (open symbols) "subbands" have 
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been folded onto the same side. From the single decade 
of data showing the increase of the density of states with 
growing distance from /x, we can infer the effective gap ex- 
ponents s e ff ~ 2.9 for / = 0.1 [Fig. ??(a)], and s e ff ~ 2.2 
for / = 0.4 [Fig. ??(b)]. The apparent dependence on 
/ arises because upon increasing the filling the flux lines 
have to accomodate more and more with the underly- 
ing random pin distribution, whereby correlations are 
destroyed. Of course, we cannot be sure that we have 
reached an asymptotic /-independent exponent with our 
simulations. It nevertheless seems clear that the mean- 



field formula (3.7) is inadequate, and that a power law 
of the form (1.3) provides a better description of the soft 
Coulomb gap, with a maximum gap exponent 



*cff X (/ « 1) 



(3.9) 



As is implied by our notation, we do not consider these 
gap exponents as universal quantities; rather, they should 
be viewed as effective exponents describing the shape of 
the Coulomb gap in some specific energy range. 

The details of the asymptotic analytical form of g(e) 
in the limit e — > /i are not our main concern. We are 
most interested in the striking shape of the density of 
states, displaying a very marked and wide correlation- 
induced "Coulomb" gap, because this has a dramatic 
qualitativ e effe ct on the current-voltage characteristics 
(see Sec. [II D| ) . One of the most important features of 
the pseudo-gap of Fig. ?? in fact rests in its remarkable 
persistence as the interaction range, set by the London 
penetration depth A, becomes finite. Figs. ??-?? depict 
a selection of results for the single-particle density of 
states for decreasing A, with various fillings / (and the 
other parameters chosen as before). As can be seen in 
Fig. ??, the situation with X/d = 5 is practically indis- 
tinguishable from the previous case of infinite interaction 
range, in overall appearance, width and general distribu- 
tion of pinning energies [cf. Fig. ??(b) and Fig. ??(b) 
for identical / = 0.4]. The overall curves ar e si mply 
shifted along the energy axis according to Eq. ( |3.8| ) (see 
also Fig. ?? below). In Fig. ?? the more extreme cases 
with (a) / = 0.1, with a very broad and flat pseudo- 
gap, and (c) / = 0.8 are shown as well. In the latter 
case, g(fj) is manifestly non-zero. This is due in part to 
the finiteness of A, thereby subjecting the system more 
and more to the underlying pinning potential random- 
ness. In addition, the fluctuations in the chemical poten- 
tial grow considerably as / is increased, thus obscuring 
the vanishing g([i) by sampling over the energies Ci [p2[ . 
Nevertheless, it is clear that the Coulomb gap closes: its 
width becomes narrower and g(fJ.) cx A~ 2 grows, as may 
be seen by comparing Figs. ??(a) and ??(a) for / = 0.1, 
Figs. ??(a) and ??(a) for / = 0.2, Figs. ??(b), ??(b), 
and ??(b) for / = 0.4, and Figs. ??(c) and ??(b) for 
/ = 0.8. Even for A ~ d < ao there remains a consider- 
able suppression of g((i) itself, accompanied by a strong 
depletion of the density of states up to |e|/2eo ~ 0.2. In 



and the effective gap exponent for X/d — 1 and / = 0.1, 
for example, is ab out s c ff « 1.2. Notice that in accord 
with our estimate ( 2.19| ), the chemical potential roughly 
scales as cx f{X/d) z as filling and interaction range are 
varied. 



C. Finite size effects 

We now discuss finite size effects, i.e.: the dependence 
on the number of defect sites Ad used in the Monte Carlo 
simulations. We have performed extensive studies for the 
cases Ad = 100, 200, 400 (for which all the previous plots 
are depicted), and 800, averaging over k — 16000/Ad 
runs with different initial defect configurations (thus sup- 
plying an equal number of data points) . The total width 
and shape of the density of states, as well as the appear- 
ance of the Coulomb gap, remained practically identical 
upon changing N-q in the above range. The sole effect 
is a chemical potential shift, or equivalently, the total 
energy per particle, for the situations with long-range 
interactions. For X/d = 1, ^(Ad) is independent of Ad 
for Ad > 200. In Fig. ??, the mean values for the ob- 
tained chemical potentials are plotted as a function of 
1/Ad for the runs with X/d = 5 and several filling frac- 
tions / < 0.5. As is to be expected, upon increasing the 
number of sites and consequently interacting particles, 
for this very long-range repulsive forces /i(Ad) grows, 
but the data approach finite values in the limit — > oo, 
close to the open symbols on the 1/Ad = axis. These 
latter points were obtained by shifting the results for 
/x(Ad = 400) of the A — > oo simulations (obtained via 
the Ewald summation procedure) according to Eq. ([3.8]), 
with X/d = 5. Thus the size dependence of the chemical 
potential is easily understood; and there seem to be few 
other finite size artifacts except that we cannot sample 
ener gies closer to fi than e w V(L±) « 2eK (L±/ X) (see 
Sec. [IB). Note that /x oc / in Fig. ??, which is equivalent 
to B ph H in the flux-line language. 



D. Current— voltage characteristics in the 
variable— range hopping regime 

We now determine how the depletion in 17(e) near the 
chemical potential affects vortex transport properties. 
An in-plane current JIB induces a Lorentz force per 
unit length fL perpendicular to J, acting on all the flux 
lines: 



<f>0 



z x J 



Accordingly, we add the term 



an intermediate range of e, one may still use Eq. (1.3) 



5F N [{r z }} = - 



/ dz^fj 

JO ,, 



Xi{z) 



(3.10) 



(3.11) 
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to our model free energy ( |2.l| ). As stated above, in the 
boson picture this additional term represents an electric 
field E = zx 31c acting on the particles carrying charge 
4>o (see Table |) . This fictitious quantity E should not 
be confused with the true electric field £(J) (the voltage 
drop induced by the current through the movement of 
flux lines). In the spirit of the thermally assisted flux- 
flow (TAFF) model of vortex transport, |53| the super- 
conducting resistivity (i.e.: the conductivity in the boson 
representation) p = £ / J may be written as 



p~ Po exp [~U B (J)/T] 



(3.12) 



where po is a characteristic flux-flow resistivity, and Ub 
represents an effective barrier height. Unlike the original 
TAFF models, Ub(J) act ually diverges as J — > in the 
Bose glass phase [see Eq. ( |l.l| )1. 

Consider first the regime where the motion of a single 
vortex is unaffected by the other flux lines in the sam- 
ple. We shall see that this is true only in an intermediate 
current regime J\ < J < J c ||. Driven by the external 
current J, a flux line will start to leave its columnar pin 
by detaching a segment of length z into the defect-free 
region, thereby forming a half-loop of transverse size r. 
The free energy price for this process consists of the elas- 
tic energy ss eir 2 /z, and the lost pinning energy ss Uqz, 
as the first vortex to move will be the most weakly bound 
one. Note that we have to re-introduce the average pin- 
ning potential Uq — ([/&) here. Upon taking into account 
the work done by the Lorentz force, the free energy of 
the loop is approximately Iq] 



SFi(r, z) sa hr 2 / z + Uqz - f L rz 



(3.13) 



Optimizing 5Fi(r, z) for /l = first, we see that for the 
saddle-point configuration 



z*^r*^h/U 
while for finite currents 



(3.14) 



r*^cU /4> J , (3.15) 
which yields the saddle-point free energy 

SF* w cl\ /2 ul /2 /4> J , (3.16) 



which we may identify as the energy barrier in Eq. ( 3.1 2| ) 
for the half-loop nucleation process || . For a sufficiently 
low current J\ , the half-loop will typically extend to the 
nearest-neighbor pin, namely when r* w d, and hence 

Ji = cUo/fod . (3.17) 

The flux line will then form a double-kink of width 

wk = dy/^/Uo , (3.18) 

which costs a free energy 2_Ek, where 



Ek = d\/iiU 



(3.19) 



These expressions allow us to cast the current-voltage 
characteristics for J\ < J < J c into the form H 



S^PoJexp [-(E K /T)(J 1 /J)] 



(3.20) 



In the half-loop regime, the intervortex repulsions thus 
do n ot affect the transport exponent p e g = 1, cf. 
Eq. (1.1). In recent experiments on BSCCO samples, 
exponents p c s ~ 1 have been found for intermediate cur- 
rents, supporting the half-loop nucleation picture |^5| , |l6| . 

For J < J\ we thus have to take the configurational 
limitations imposed by the other vortices into account. 
The most important thermally activated excitation will 
now be a double superkink (Fig. ??): the flux line will 
send out a tongue-like pair of kinks to another possi- 
bly very distant defect, which has a favorable pinning 
energy, thus optimizing the tunneling probability. This 
is the flux line analogue of variable-range charge trans- 
port in disordered semiconductors |Pfl|. The cost in 
free energy for such a configuration of transverse size 
R and extension Z along the magnetic-field direction 
will consist of two terms: (i) the double-superkink en- 
ergy 2Ek(R) — 2Eyi(d)R/ d stemming from the elastic 
term, and (ii) the difference in pinning energies of the 
highest -energy occupied site, ~ p, and the empty site 
at distance R with ej = fi+ A(i2). Thus the free energy 
difference with respect to the situation without kinks and 
external current is 



6F SK » 2E K R/d + ZA(R) - f L RZ , 



(3.21) 



which generalizes Eq. (4.13) of Ref. || for the case of 
non-interacting flux lines. The concentration available 
states as a function of R with D dimensions trans- 
verse to B (here D = 2) on the one hand equals 

d D J/T +A ^ 5( e )^ e i an d on the other hand is simply given 
by (d/R) D ; thus A(i?) is to be determined from the 
equation 



/ 



M-A(-R) 



g(e)de = R 



-D 



(3.22) 



Optimizing now first for vanishing current J = gives 
the longitudinal extent Z* of the superkink as a function 
of its transverse size R*. 



2E K /d 
(dA/dR)! 



(3.23) 



Upon balancing the last term in Eq. (3.21) against the 
optimized sum of the first two, one arrives at 



J^ /c« A(R*)/R* 



(3.24) 



which through inversion yields a typical hopping range 
R*(J). Inserting back into (3.21) finally yields the result 
for the optimized free energy barrier for a jump, 
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8F£ K (J) « (2E K /d)R*(J) 



(3.25) 



which we ident ify w ith the current-dependent activation 
energy in Eq. ( 3.12 ) 



£ sa poJ"exp [-(2£^/Td)i?*(J)] 



(3.26) 



Consider now a power-law form for the distribution of 
pinning energies, g(e) = n\e — /i| s . Then from Eq. (3.22) 



A(R) 



V(«+i) 



and ( |3.24| ) yields 
R*(J) = 



R -D/(s+l) 



. , \ l/(Z3+s+l) , x p 

s + 1 \ / c ■ 



(3.27) 



where 



V - 



s + 1 



D + s + 1 
This immediately leads to 

<SF* K (J) = 2E K (Jo/J) p 



(3.29) 



(3.30) 



which is of the form (1.1) with the transport exponent 

(3.31) 



(3.29) and the current scale 

J « c/0 o « 1/(s+1) rf 1/p 



For a constant density of states, s = and k = g(/i), 
one recovers the non-interacting results, namely the 
Mott variable-range hopping exp onen t in D dimensions 
Po = 1/(15 + 1), and expression (|l^) for J g. Using 
the mean -field r esult for the gap exponent, sq = D/a—l 
[Eq. Q)], Eq. ( ^29|) gives p = 1/(1 + 0"), i.e.: p^ 1 for 
purely logar ithmic interactions. In the light of our results 
in Sec. IIIB| , however, we have to expect that p e g will ac- 
tually be smaller; namely with Eqs. ( [3.9] ) and ( 3.29| ) we 
find 



max 
PcS 



2/3 



(3.32) 



in two dimensions, for A —> o o and small fillings. 

The numerical result ( 3.32] ) for the variable-range hop- 
ping transport exponent of logarithmically interacting 
particles in two dimensions is consistent with the anal- 
ysis by Fisher, Tokuyasu, and Young, |^5| although the 
underlying theoretical arguments appear to be different. 
One might also question whether the "gauge glass" model 
studied in Ref . [p5| , with its quenched random vector po- 
tential, correctly captures the physics of a finite density 
of, say, positive vortices subject to a random pinning po- 
tential. 

Using the general relations (3.22)-( [3.26| ) and the 
Monte Carlo results of Sec. Ill B for the single-particle 
density of states, we can numerically evaluate the 



current-voltage characteristics for any form of g(e). Re- 
sults derived from the distributions of interacting pin- 
ning energies in Figs. ??-?? are depicted in Fig. ?? 
as double-logarithmic plots of the effective hop size 
R*(J) oc SFg K (J) = U b(J) as function of the current J; 



according to Eq. ( 3.26 ), the slope immediat ely y ields the 
effective transport exponent p e g [cf. Eq. ( [3.28] )]. Each 
plot corresponds to a certain filling fraction /, and sum- 
marizes the results for various values of the parameter 
X/d. In Fig. ??(a) for / = 0.1 the curve corresponding 
to long-range repulsions (X/d — 5, which turns out to 
be practically indistinguishable from the case A — > oo) 
yields p c g w 0.70 ± 0.05, in accord with ( 3.32j ). Upon 
decreasing the interaction range, p R becomes smaller, 
Pcff ~ 0.55 ± 0.05 for X/d = 2 and p cS « 0.50 ± 0.05 for 
X/d = 1, and finally reaches the non-interacting expo- 
nent po = 1/3 for A — > 0. Note that there appears as 
well a considerable shift in the prefactor for the power 
law which also enhances the effectiveness of the vor- 
tex pinning, when the repulsive interactions are turned 
on. Long-range repulsive interactions thus reduce vor- 
tex voltages oc exp[— 2E-&R* (J) / dT] by several orders of 
magnitudes with respect to the non-interacting situation. 
E.g., for log j w 1.5 the exponent in ( [3.26] ) is about ten 
times smaller for X/d = 5 as compared to A — > 0, and 
hence the resistivity is reduced by a factor of w 10~ 5 
[note that \i as function of A does actually not vary very 
much, as e oc A~ 2 , see Eq. ( 2 . 1 9|) 1 - 

For higher fillings /, this collective pinning mechanism 
becomes relatively weaker, as the flux lines increasingly 
have to accomodate the underlying random pin distribu- 
tion, and correlations are destroyed (see the discussion in 
Sec. IIIB). This can be seen in Fig. ??(b) and (c); for 



/ = 0.8 all the curves with different finite A are essen- 
tially parallel to the the one for A — > 0, i.e.: p e ff ~ 1/3- 
However, there is still a considerable offset due to the 
change in the prefactor of the po wer l aw, and a reduc- 
tion of the exponential factor in ( |3.26 ) by about three, 
and hence p still about twenty times smaller than for 
the case without interactions. We remark that of course 
the data for the very lowest currents in Fig. ?? are less 
reliable, for the finite size of our system becomes impor- 
tant when R*(J) w L±, which for Nd = 400 happens for 
log[#*( J)/d] w 1.3. 

In at least one recent experiment on BSCCO, reported 
by Konczykowski, Chikumoto, Vinokur, and Feigel'man, 
|15| an effective transport exponent p c g in the variable- 
range hopping regime clearly different from the non- 
interacting po — 1/3 was seen. The relevant parame- 
ters of the measurement at T w 60K and H = 300Oe 
depicted for B = 0.2T in Fig. 4 of Ref. |]| are / = 
B/B^ « 0.15, and, using A(0) w 1400A, X(T)/d w 1.6. 
Across about a half-decade the experimental current- 
voltage characteristics could be described by an effective 
exponent p c g w 0.57, which agrees very well with our 
result p c s ~ 0.55 for / = 0.1 and X/d = 2 (note that 
both sets of parameters correspond to the same ratio 
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A/ao ~ 0.6). For the experiment with much higher irradi- 
ation dose, B,p = 2T, corresponding to much lower fillings 
/ w 0.015 and A/d w 5, however, the variable-range hop- 
ping exponent turned out to be lower, namely p c g « 0.39, 
p6| contrary to the expectations from our above calcula- 
tions. A possible explanation is that for such low values 
of /, see also Ref. jl4|] where p c s ~ 1/2 was measured for 
/ w 0.0065, one has to reconsider the above analysis and 
take very rare fluctuations in the spatial distributions of 
rods into account, leading to p c ff = 1/2, at least for the 
non-interacting case j8|]. 

Finally, we briefly address the case of samples with fi- 
nite thickness L. At the very lowest currents J — * 0, 
the current-voltage characteristics will event ually be- 
come Ohmic, § namely when L rs Z*, cf. Eq. ( 3.23 ); for 
the power-law density of states ( [L.3| ) studied above the 
typical longitudinal extent of a supcrkink is 



Z*(J) = {2E K /d)(c/<f> J) 



(3.33) 



and hence the crossover current scale to Ohmic resistance 
becomes 



J L = (2E K /d){c/<f> Q L) . 
Substituting this into Eq. fl3.30| ), we find 

6F£ K (L)/T=(L/L y 
with the characteristic length scale 

L m n 1/{s+1) (d/2E K ) D T 1/p . 
The ensuing resistivity law 

p = Po cxp [-(L/L ) p ] 



(3.34) 
(3.35) 
(3.36) 
(3.37) 



with its unusual thickness-dependence is the direct ana- 
logue of the variable-range hopping conductivity formula 
for disordered semiconductors, |Q generalizing Mott's 
law to the case of long-range interactions. The latter is 
recovered when s — > and hence in D = 2 dimensions 
p = po = 1/3 and L = g(p)(d/2E K ) 2 T 3 §. 



IV. SUMMARY AND DISCUSSION 

We have studied spatial correlations and the inter- 
acting density of states of flux lines pinned to colum- 
nar defects in the Bose glass phase ||. Our procedure 
was to map the vortex problem onto an equivalent two- 
dimensional disordered boson system at zero tempera- 
ture, and examine the latter using an established 
Monte Carlo simulation algorithm ]l7| , p2| -p4]| . We ex- 
pect our simplified Hamiltonian to apply for B < B* 



[Eq. d2.ll] )] and "low" temperatures T < T x [Eq. (gT2J)], 
which may in fact extend throughout a large fraction 
of the Bose glass phase, up to temperatures very near 
to the irreversibility line. The resulting distribution of 



pinning energies g(e), which explicitly takes the possi- 
bly long-range interactions into account, was then used 
to calculate the vortex transport characteristics in the 
varia ble-r ange hopping regime at low currents J < J\ 
[Eq. ( 3.17 )], by generalizing the optimization procedure 
of Ref~l§] . 

For any A > ao and low fillings / < 0.3, we have found 
that the flux lines form a highly correlated structure, the 
corresponding S(q) displaying a marked, but broad peak 
at qao ~ 2n. This spatial configuration may be under- 
stood as a highly distorted triangular lattice, with a large 
number of topological defects in the vortex coordination 
shells imposed by the underlying random distribution of 
columnar pinning sites. The peak in S(q) is destroyed for 
larger fillings, and/or for weaker interactions A < ao. As 
for the Coulomb glass, |2^] a peculiar spatial clustering 
of those sites with pinning energies near the chemical po- 
tential fi was found, the empty low energy sites and filled 
high energy sites occupying roughly complementary re- 
gions of space. As opposed to the case of localized charge 
carriers in disordered semiconductors, fllTj the above spa- 
tial structures can be directly investigated in experiment, 
and the first such measurements have already been re- 
ported 12 13]. Quite unambiguously, highly correlated 



vortex distributions were found, and thus the flux line re- 
pulsions need to be accounted for |D|]. Moreover, we have 
demonstrated how measured data on columnar pin and 
vortex positions may be utilized to actually predict the 
corresponding density of states and estimate the current- 
voltage characteristics Q . 

In the distribution of pinning energies (single-particle 
density of states), for A > ao the interactions lead to 
the formation of a remarkably prominent and persisting 
"Coulomb" gap near the chemical potential /i separating 
the occupied and empty states. This soft pseudo-gap can 
at least approximately be described by a power law ([L"s|), 
characterizing the depletion of states upon approaching 
\x. With increasing filling / for fixed interaction range A, 
the gap exponent becomes smaller and the pseudo-gap 
eventually closes as the vortices have to accomodate with 
the underlying randomness, thus losing their spatial cor- 
relations. With our simulations, we found effective gap 
exponents s c g 1 for A » d up to values s e ff ~ 3 for 
A — > oo and /< 1, which demonstrates that correlation 
effects are in fact much stronger than predicted by the 



qualitative mean-field estimate (3.7) 



These equilibrium results were used as a basis to dis- 
cuss both the half-loop nucleation region at intermediate 
currents, as well as the variable-range hopping transport 
regime at very low currents, where the vortices move by 
forming double-superkinks to favorable, possibly distant 
sites |8|]. For the former case, we found that p e g = 1 as 
in the non-interacting situation. In the variable-range 
hopping regime, however, the depletion of low-energy 
states according to Eq. (1.3) implies a change of the ef- 
fective transport exponent p c g to higher values, render- 
ing the collective pinning of vortices to columnar defects 
much more effective when A > ao- For A = d we have 
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found p e g w 1/2, while for very long-range repulsion 
A — > oo values up to p e ff ~ 2/3 may be reached, similar 
to the scaling analysis by Fisher, Tokuyasu, and Young 
p5| . In at least one measurement an effective exponent 
p cf f w 0.57 was found for parameter values well in accord 
with our simulations [fj5"| . 

Therefore, we may speculate that the remarkable cor- 
relation effects reported in this paper, which are induced 
by the long-range intervortex repulsions, could in fact 
be utilized for designing future superconducting materi- 
als in order to reduce dissipative flux transport as much 
as possible. For example, one could specifically "tailor" 
samples to obtain large values of A/ao- The remarkably 
strong pinning predicted for splayed columnar defects 
fl9f might be even further enhanced by interactions. 
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FIG. ??. Columnar defects and pinned flux lines, with 
filling fraction / < 1. 

FIG. ??. Schematic phase diagram of flux lines inter- 
acting with columnar defects that are aligned along the 
magnetic field direction. 

FIG. ??. Double-superkink configuration, which is the 
most important excitation relevant for the variable-range 
hopping transport regime of vortices. 

FIG. ??. Schematic current-voltage characteristics in 
the Bose glass phase. 

FIG. ??. Spatial positions of unoccupied columnar 
pins (open circles) and columnar pins occupied by mag- 
netic flux lines (filled circles), as obtained from typical 
simulations (with randomly chosen initial vortex config- 
uration) using an identical set of underlying defect coor- 
dinates (iVo = 400). The parameters are: w/2e = 0.1, 
A — > oo (i.e.: logarithmic interactions), and (a) / = 0.1, 

(b) / = 0.2, and (c) / = 0.4. 

FIG. ??. Vortex structure function S(q), obtained by 
averaging over the initial and final states of a series of 
k = 40 different defect configurations with A D = 400, 
w/2e = 0.1, X/d = 5, and (a) / = 0.1, (b) / = 0.2. S(q) 
for the initial flux line distributions is the long-dashed 
line for randomly chosen positions, and the dashed line 
for the distorted triangular lattice; the final vortex dis- 
tributions are depicted as the dot-dashed line for the 
configurations stemming from random initial states, and 
the solid line originating in the distorted triangular lat- 
tice positions. Notice the shift in the peak position of the 
dashed curves in (a) and (b). 

FIG. ??. Dependence of the vortex structure func- 
tion S(q) on the filling /. The results were obtained 
by averaging over k = 40 different defect configurations 
(with TVd = 400 and random initial occupations), where 
w/2eo = 0.1 and X/d = 2; solid: / = 0.1, dot-dashed: 
/ = 0.2, dashed: / = 0.4, and long-dashed: / = 0.8. A 
broad peak at qao w 2n is found only for / < 0.4. 

FIG. ??. Voronoi analyses (Delaunay triangulations) 
of the vortex positions shown in Fig. 5; i.e.: Nu = 400, 
w/2e = 0.1, A ^ oo, and (a) / = 0.1, (b) / = 0.2, 

(c) / = 0.4. The vortices with sixfold coordination arc 
depicted as full circles, while topological defects with co- 
ordination number N c — 4 are shown as open circles, 
N c = 5 as open diamonds, N c = 7 as encircled asterisks, 
and N c = 8 as encircled crosses. The histograms depict 
the relative abundancy h(N c ) for the occurence of the 
coordination numbers iV c = 4, 5, 6, 7, 8, respectively. 

FIG. ??. Spatial clustering of both occupied (filled cir- 
cles) and empty (open circles) sites with energies near the 
chemical potential, for the final configuration of Fig. 5(c), 
i.e.: N D = 400, / = 0.4, w/2e = 0.1, and A -> oo. The 
small symbols depict the complete distribution of vortices 
and pins, while the large circles correspond to those sites 
whose energies ej lie in the range (a) \a — fi\/2e < 0.6, 
and (b) |ej — yu|/2e < 1.2, respectively. 



FIG. ??. Normalized distribution of pinning energies 
G(E') — 2ead 2 g(e) as function of the single-particle en- 
ergies E 1 = e'/2eo, averaged over k = 2 x 40 runs with 
N D = 400, w/2e = 0.1, A ^ oo, and (a) / = 0.2, (b) 
/ = 0.4; /x'/2eo ~ —1.3 (marked by the arrow) in both 
cases. Open circles: results of the simulations starting 
from random initial configurations, filled triangles: site- 
energies from runs where a distorted triangular lattice 
was used as initial state. 

FIG. ??. Double-logarithmic plots of the normal- 
ized distribution of pinning energies G(E) = 2eod 2 g(e)vs. 
\E — 2? M | = |e— /x|/2eo, averaged over k = 2x40 runs with 
N D = 400, w/2e = 0.1, A ^ oo, and (a) / = 0.1, (b) 
/ = 0.4. Circles and triangles refer to the choices of ran- 
dom positions or a distorted triangular lattice as initial 
configurations, as in the previous figure. Filled and open 
symbols refer to single-particle energies in the occupied 
(n, = 1) and empty (n, = 0) pseudo-bands, respectively. 

FIG. ??. Normalized distribution of pinning energies 
G(E) = d 2 g(e) as function of the single-particle energies 
E = e/2e , averaged over k = 2 x 40 runs with iV D = 400, 
w/2e = 0.1, X/d = 5, and (a) / = 0.1: ^/2e ~ 10.9, 
(b) / = 0.4: ^/2e w 46.9, (c) / = 0.8: ^/2e w 96.8. 
The symbols have the same meaning as in Fig. 10. 

FIG. ??. As in Fig. 12, but for X/d = 2 and (a) / = 0.1: 
^/2e « 1-74, (b) / = 0.4: ^/2e w 9.07. 

FIG. ??. As in Fig. 12, but for X/d = 1 and (a) / = 0.2: 
fi/2e w 0.72, (b) / = 0.8: ^/2e w 5.33. 

FIG. ??. Chemical potential = /x/2eo as function of 
the system size (7V D = 100, 200, 400, 800), averaged over 
k = 16000/A^d runs (with random initial vortex distribu- 
tions), for u>/2eo = 0.1 and X/d = 5; squares: / = 0.1, 
triangles: / = 0.2, circles: / = 0.3, inverted triangles: 
/ = 0.4, diamonds: / = 0.5. The open symbols plot- 
ted at 1/Nj) = represent the corresponding results for 
ft = fj,'(Nu = 400) + 27r/(A/d) 2 obtained with the chem- 
ical potential results /x' for the purely logarithmic poten- 
tial (A — > oo). 

FIG. ??. Double-logarithmic plots of the exponential 
factor R*(J)/d oc 5F£ K (J) for variable-range hopping 
vs. the reduced current j = J(j} d/2eoc, derived from 
the numerically obtained density of states, averaged over 
k = 40 runs with Ad = 400, random initial configura- 
tions, and w/2e = 0.1, for different interaction ranges: 
diamonds: A — > 0, circles: X/d — 1, squares: X/d = 2, 
triangles: X/d = 5; (a) / = 0.1, (b) / = 0.4, and (c) 
/ = 0.8. 
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